Effect of gold nanoparticles distribution radius on photothermal therapy efficacy

Lasers are used in various fields, however, in the medical field, they are mainly used for incision or chemotherapy. Photothermal therapy (PTT) is an anti-cancer treatment technique that uses lasers and the photothermal effect to increase the temperature of tumor tissue and induce its death. In this study, the therapeutic effect of PTT using gold nanoparticles as a photothermal converter was analyzed numerically for the occurrence of squamous cell carcinoma inside a skin section consisting four layers. Numerical modeling was implemented to calculate the temperature distribution inside the biological tissue while varying the distribution radius of gold nanoparticles in the tumor tissue, the number of injections, and the intensity of the irradiating laser. For the given situation, the optimal treatment effect was observed when the distribution radius ratio of the injected gold nanoparticles (GNPs) was 1, the number of injections was 7, and the intensity of the irradiated laser was 52 mW. Three apoptotic variables were used to quantitively evaluate the effect of PTT in each case and thus suggest the optimal treatment effect. However, although the temperature range at which apoptosis occurs is known, the maintenance of that temperature range is still under research and the temporal influence of apoptosis remains to be determined.


Effect of gold nanoparticles distribution radius on photothermal therapy efficacy Donghyuk Kim , Jeeyong Paik & Hyunjung Kim *
Lasers are used in various fields, however, in the medical field, they are mainly used for incision or chemotherapy. Photothermal therapy (PTT) is an anti-cancer treatment technique that uses lasers and the photothermal effect to increase the temperature of tumor tissue and induce its death. In this study, the therapeutic effect of PTT using gold nanoparticles as a photothermal converter was analyzed numerically for the occurrence of squamous cell carcinoma inside a skin section consisting four layers. Numerical modeling was implemented to calculate the temperature distribution inside the biological tissue while varying the distribution radius of gold nanoparticles in the tumor tissue, the number of injections, and the intensity of the irradiating laser. For the given situation, the optimal treatment effect was observed when the distribution radius ratio of the injected gold nanoparticles (GNPs) was 1, the number of injections was 7, and the intensity of the irradiated laser was 52 mW. Three apoptotic variables were used to quantitively evaluate the effect of PTT in each case and thus suggest the optimal treatment effect. However, although the temperature range at which apoptosis occurs is known, the maintenance of that temperature range is still under research and the temporal influence of apoptosis remains to be determined. www.nature.com/scientificreports/ different breast cancer models. PTT was performed in two cancer models using 50 μg/ml of non-toxic MXene, and a significant decrease in tumor survival was observed after treatment along with a significant increase in reactive oxygen species. This confirms that PTT with single and multilayer MXene can regulate tumor growth by expressing apoptosis through temperature elevation. Wang et al. 30 identified the optimal temperature distribution of the tissue injected with gold nanoparticles (GNPs) to achieve better results in clinical applications. The heat generation in the biological tissue irradiated with the laser was calculated using the Monte Carlo method, and the temperature distribution inside the tissue was calculated using the two-energy equation. In addition, the Arrhenius equation was used to determine the permanent thermal damage of tissues. A parametric study was conducted by varying the intensity of the laser, the volume fraction of GNPs, the heating and cooling time, and the irradiation position, and it was found that controlling the heating and cooling time could effectively prevent overheating of the skin surface and thermal damage to internal tissues. In addition, it was confirmed that lower volume fraction of GNPs achieved better PTT effects when using the ring heating strategy. As a result, PTT requires a combination of biological, optical, nano, and heat transfer research, and there is still significant scope for improvements. In the field of heat transfer, theoretical research has focused on evaluating the temperature distribution in biological tissue. However, this does not reflect the degree of apoptosis that is commonly used in biology, and simply checks for thermal damage to the tissue through the Arrhenius thermal damage model. In addition, in most of the previous studies, numerical analysis was performed after assuming that the GNPs were uniformly distributed. Furthermore, in the field of medicine and biology, phenomenological studies have only confirmed results under certain predefined conditions and have not determined the optimal conditions for PTT.

List of symbols
Therefore, this study investigated the conditions under which the PTT effect is maximized by changing the distribution radius of GNPs in the tumor tissue. The behavior of the laser particles in the medium was analyzed using the Monte Carlo method and the optical properties of the GNPs were calculated using the Discrete Dipole Approximation method. In addition, the temperature distribution in the medium was calculated based on the heat diffusion equation for a simulated PTT situation. Finally, the effectiveness of PTT was confirmed through the apoptotic variable, which quantitatively judges the degree of maintenance of the apoptosis temperature range during the treatment time and the amount of thermal damage to the surrounding normal tissue, as proposed by Kim et al. 31 and Kim & Kim 32 , and the optimal treatment conditions were suggested.

Governing equation and numerical process propagation
Calculation of light absorption and heat distribution of medium. In this study, the Monte Carlo method was applied to simulate laser irradiation of biological tissue. The Monte Carlo method is a technique for analyzing the behavior of laser particles in a medium whereby the absorption and scattering behavior are analyzed simultaneously to calculate the final light absorption distribution 33 . By calculating the travel distance S, deflection angle θ, and azimuth angle ψ at each time step using random numbers ξ, the path of travel in the medium for one particle can be determined. The final distribution of light absorption in the medium can be determined through probability distribution analysis by repeatedly calculating the path for each of the number of particles set initially. The relevant equations are shown below.
Equation (1) is the formula for calculating the distance traveled by a given particle in a single time step based on the total attenuation coefficient of the medium (μ ext ) and a random number 33 . The total attenuation coefficient can be expressed as the sum of the absorption coefficient (μ abs ) and reduced scattering coefficient (μ ' sca ). Equations (2) and (3) can be used to calculate θ and ψ, respectively 33 . For θ, there are two applicable different formulae, depending on the value of the anisotropy factor (g), which is a dimensionless number indicating the amount of energy retained in the forward direction after a single scattering. Both θ and ψ are also determined using a randomly selected random number.
In the case of a laser beam, the amount of laser heat absorbed by the medium depends on the profile being irradiated. The profile of the laser is broadly divided into gaussian and top-hat, which can be irradiated into the medium in the form of Eqs. (4) and (5), respectively 33 .  www.nature.com/scientificreports/ Once the distance and angle moved by the particle are calculated, the percentage decrease in energy W per time step can be calculated based on the optical properties of the medium, as shown in Eq. (6) 33 . The particle continues to move until, finally, its energy converges to 0. Equations (1-6) are repeatedly calculated for each of the set number of particles. Following that, the final light absorption distribution in the medium can be calculated by calculating the probabilistic particle distribution as shown in Eq. (7), where φ z represents the energy density of the photon in the depth direction and φ rz is the absorbed photon probability density function 33 . In addition, N denotes the number of initially selected photons, and i r and i z are the indices for the grid elements in the r and z directions, respectively.
Finally, the light absorption distribution in the medium (Eq. (7)), obtained using the Monte Carlo method, is used in calculation of the temperature distribution in the medium based on the heat diffusion equation, as shown in Eq. (8), and the finite element technique (Eq. (9)) 34 .
Optical properties of gold nanoparticles and biological tissue. The absorption and scattering coefficients in the Monte Carlo method require that the optical properties of nanoparticles and medium must first be calculated. Accordingly, in this study, the discrete dipole approximation (DDA) method was used to calculate the optical properties of GNPs, which were the photothermal converters in this study, and the biological tissues into which the GNPs were injected 35 . This method enables calculation of the optical properties of nanoparticles of various shapes, whereas Mie theory is limited to conventional spherical or elliptical nanoparticles 36 . The DDA method used in this study is detailed below.
First, the polarization vector P, which represents the dipole moment in the unit volume, is calculated as shown in Eq. (10) from the polarizability α and the local electric field E 37 . The variable E can be calculated using Eqs. (11) and (12), where r and k represent the position vector and wavenumber, respectively, and A is the interaction matrix between the dipoles, calculated as in Eq. (13) 37 . If i and j are equal, the interaction matrix A in Eq. (13) can be simplified to α i −1 .
The calculated P can then be used to calculate the total attenuation C ext , absorption C abs , and scattering cross sections C sca , as shown in Eqs. (14)(15)(16), where * indicates the complex conjugate 37 . www.nature.com/scientificreports/ This allows the final absorption Q abs , attenuation Q ext , and scattering efficiencies Q sca of the medium to be calculated as shown in Eq. (17), where r eff is the effective radius of the nanoparticle (from Eq. (18)) and V is its volume 35 .
Once the optical efficiency of a nanoparticle is calculated, the optical properties of the medium, considering the presence of the nanoparticles, can be calculated using the relation proposed by Dombrovsky et al. 38 . The relation is determined by the volume fraction of nanoparticles in the medium f v , their respective optical efficiency, and r eff , and is calculated as in Eqs. (19,20) 38 . Following the calculation of optical properties of the nanoparticles, the final optical properties of the medium into which the nanoparticles are injected can be calculated as the sum of the calculated optical properties of the nanoparticles and the optical properties of the medium itself, as shown in Eq.(21) 38 . This makes it possible to calculate the optical properties of the entire medium in the presence of nanoparticles, which can then be applied to the Monte Carlo method.
Verification and comparison of numerical model. To validate the numerical model presented in this study, a comparative analysis was performed with the results presented by Gnyawali et al. 39 . The study assumed the presence of a spherical tumor with a radius of 5 mm inside a cylindrical biological tissue with a radius of 30 mm and a depth of 30 mm, as shown in Fig. 1a. A laser with a top-hat profile, a wavelength of 805 nm, and a beam radius of 15 mm was utilized, irradiating the top surface with an intensity of 1 W/cm 2 for 600 s. The initial temperature of the tumor and biological tissues was assumed to be 0 °C. The analysis was performed using the numerical modeling presented in this study under the same conditions as in 39 .
The Fig. 1b shows the results of the numerical simulation verification. The results of 39 are shown as a red line, and the results of the numerical modeling presented in this study are shown as black scatter points. A root mean square error of 0.6775 was obtained compared to the results of previous studies and, as shown in the graph, the numerical analysis model proposed in this study fits well with previous studies. In addition, a biomimetic phantom was used to validate the numerical model and further demonstrated its validity 32 . www.nature.com/scientificreports/ Numerical geometry and calculation conditions. A situation in which squamous cell carcinoma (SCC) occurred inside human skin composed of four layers was numerically modeled, as shown in Fig. 2.
Assuming that SCC with a diameter of 4 mm and a length of 2 mm occurred in the center of a skin layer with a diameter of 20 mm and a depth of 10 mm, a continuous wave laser with a Gaussian distribution, a wavelength of 1064 nm, and a diameter of 4 mm was irradiated in the perpendicular direction to the skin surface. The laser irradiation time was fixed at 600 s. The length and thermal/optical properties of each skin layer and tumor are shown in Table 1 [40][41][42][43][44][45][46] .
In this study, the temperature distribution of the tumor and surrounding normal tissues was investigated while changing the distribution radius of GNPs within the tumor tissue. Figure 3 shows the change in the distribution radius of GNPs in the tumor when viewed in the vertical direction.
Two situations, a single GNPs injection and GNPs deposition over seven injections, are shown as examples. As evident in the figure, the distribution radius of each injection was also varied. Equation (22) shows the distribution radius ratio of GNPs (φ drr ), defined as the ratio of the distribution radius of the GNPs in the base case ("Basis") to the current distribution radius of GNPs.
The total volume of injected GNPs was kept constant between each number of injections and the distribution radius of GNPs was decreased by 20% from the basis over five steps. For the basis, the initial radius of the injected GNPs was 1 mm for a distribution radius ratio of 1 in a single injection; the radius of the injected GNPs decreases as the distribution radius ratio decreases and the number of injections increases. The final numerical solution conditions are summarized in Table 2.
The distribution radius ratio of GNPs φ drr was varied from 0.2-1.0 at intervals of 0.2. The number of injections was increased from 1 to 7, at intervals of 1. The distribution radius was changed according to the number of injections. Therefore, there were a total of 35 different cases for the distribution radius of GNPs. Since the total amount of injected GNPs was kept constant between the number of injections, the volume fraction of GNPs in the medium due to each injection decreases as φ drr decreases and the number of injections increases.  www.nature.com/scientificreports/ The volume fraction of GNPs in the medium was defined as the ratio of the volume of the radius in which the GNPs were distributed to the volume of the entire tumor tissue, from which the optical properties of the space in which the GNPs were distributed were calculated. In addition, numerical analysis was performed by varying the intensity of the irradiated laser from 0-100 mW at 2 mW intervals. Accordingly, the temperature distribution was calculated for 1785 individual cases. Using these results, the PTT conditions that showed the optimal treatment effect were identified.

Results and discussion
Result of temperature field. To calculate the apoptotic variable proposed by Kim & Kim 32 to quantitatively indicate the effectiveness of PTT, it is first necessary to determine the temperature distribution inside the biological tissue. As mentioned above, in this study, the temperature change resulting from varying the distribution radius of GNPs inside the skin cancer tumor was evaluated. The Fig. 4 shows the temperature distribution in the YZ plane of normal and tumor tissue at 200 and 600 s after laser irradiation, respectively, for φ drr = 1, number of injections = 3, and P l = 50 mW. The laser was set to a wavelength of 1064 nm with a gaussian profile. The area inside the black box represents the tumor tissue. After 600 s had elapsed, the temperature was higher than at the time of 200 s has elapsed. Tumor tissues were confirmed to have entered the apoptosis temperature range for approximately 96% of the tumor volume at 200 s and at 94% at 600 s. However, after 600 s, the normal tissue underlying the tumor tissue reached a temperature of around 48.5 °C, which may cause thermal damage to the normal tissue.
The Fig. 5 shows the apoptosis and necrosis areas in the tissue at , depth = 1 mm, number of injections = 7, P l = 100 mW, and 150 s after the start of treatment. In the figure, the region marked in green represents the apoptosis temperature range, red represents necrosis, and blue represents the region below 43 °C. As shown in the figure, even under the same conditions, the area that satisfies the apoptosis temperature range within the tumor tissue is different depending on φ drr . For φ drr of 0.4, the percentage of satisfying apoptosis temperature range in tumor tissue was derived to be about 87%, and for φ drr of 0.8, it was derived to be about 19%. In addition, in the case of normal tissue around the tumor, most of the area is in the apoptosis temperature range, but apoptosis in normal tissue can be considered as thermal damage, so it is necessary to minimize thermal damage by adjusting the appropriate laser intensity.  www.nature.com/scientificreports/ As illustrated above, this study successfully obtained the temperature distribution inside the medium at each time step in all the considered cases using the proposed method, and quantitatively confirmed the effect of PTT by calculating the apoptotic variable (described in the next section).

Degree of apoptosis maintenance in the tumor tissue.
Using the temperature distribution in the medium calculated in Section "Result of temperature field", the apoptosis retention ratio ( θ * A ) was derived to determine the degree of retention of the apoptosis occurrence temperature range inside the tumor tissue. The variable θ * A represents quantitatively how much of the tumor tissue was maintained within the temperature range of 43-50 °C (which, as described earlier, is where apoptosis is known to occur). It is calculated as the ratio of the volume of tumor tissue corresponding to the apoptosis temperature range to the total volume of the tumor tissue (Eq.(23)) 32 . Through this variable, it is possible to quantitatively confirm how well apoptotic temperature is maintained during the treatment period.
The Fig. 6 plots θ * A as a function of φ drr and P l for the number of GNPs injections. The laser was set to a wavelength of 1064 nm with a gaussian profile. First, it was observed that P l decreased as φ drr increased for all injections, with θ * A exhibiting a clear maximum in each situation. This is because, as φ drr increases, the distribution radius of GNPs in the tumor increases, and thus the light absorption area increases, enabling the absorption of more heat and causing a higher temperature rise. Furthermore, a single GNPs injection resulted in a relatively low θ * A , which is attributed to the GNPs being concentrated in the central part of the tumor tissue and thus not enabling heating of the entire tumor area. When reviewing all cases, it was found that θ * A was maximized when φ drr = 0.8, the number of injections = 7, and P l = 68 mW, indicating that the intratumoral apoptotic temperature Total tumor volume dτ  www.nature.com/scientificreports/ is maximized in this case. Compared to the author's previous study, which was analyzed under the assumption that GNPs were uniformly distributed throughout the tumor 47 , it was confirmed that P l with the maximum θ * A was derived small. In addition, it was confirmed that the optimal value of θ * A itself was derived about 10% higher than the previous result. Compared to the situation in which GNPs are assumed to be distributed throughout the tumor tissue, a relatively small amount of GNPs is injected in the case where GNPs are partially injected, so that the range of absorbing laser heat inside the medium is reduced. Compared to the situation where GNPs are distributed throughout the tumor, the situation where GNPs are partially injected results in a relatively smaller amount of GNPs being injected, leaving less area within the medium to absorb the laser heat. This results in fewer heat-generating zones, which means less unnecessary temperature rise and more efficient temperature control inside the tumor tissue.
However, even if the temperature at which apoptosis occurs in the tumor tissue is maintained for as long as possible, this does not necessarily indicate an optimal treatment condition. the thermal damage to the surrounding normal tissues may also increase further. Hence, evaluation should consider the final treatment effect, including the thermal hazard retention value (as described later).
Thermal damage to surrounding normal tissues. To reiterate, in PTT, tumor tissue is irradiated with a laser, causing a temperature increase that kills the tumor tissue. The heating itself is concentrated in the tumor tissue, but the temperature of the surrounding normal tissue is also raised through heat conduction between it and the tumor tissue. Excessive temperature increment cause thermal damage to surrounding normal tissue, which needs to be analyzed quantitatively. Therefore, in this study, the thermal hazard retention value ( θ * H ) was used to quantify the amount of thermal damage. The variable θ * H is defined as the ratio of the volume of tissue at a temperature indicating thermal damage to the weighted sum of the volume of the surrounding normal tissue; weights are assigned to the various phenomena (e.g., protein denaturation, carbonization) expressed at the various temperature ranges in the biological tissue (Eq.(24)) 32 . The surrounding normal tissue was selected to be 50% of the length of the tumor, so the thermal damage in that area can be identified. www.nature.com/scientificreports/ The Fig. 7 plots θ * H as a function of φ drr and P l for the number of GNPs injections. The laser was set to a wavelength of 1064 nm with a gaussian profile. As expected, as the P l increases, the values of θ * H increase. It was also observed that, for the same P l , θ * H increased with increasing φ drr , which is explained by the increase in light absorption in the tumor tissue, allowing more heat to be absorbed. Furthermore, it was found that the amount of thermal damage to the surrounding normal tissue was similar regardless of the number of injections. Compared to the situation where GNPs are assumed to be uniformly distributed, it was found that less thermal damage to the surrounding normal tissue was obtained 47 . This is because there are fewer heat-generating areas within the tumor tissue, as mentioned in Section "Thermal damage to surrounding normal tissues", so less heat is conducted to the surrounding normal tissue.
PTT effect analysis and optimal treatment conditions. Finally, to analyze the effectiveness of PTT, it is necessary to consider both the degree of apoptosis maintenance and the degree of thermal damage to surrounding tissue simultaneously. Since the previously calculated θ * A and θ * H confirm the thermal effects of only tumor tissue and surrounding normal tissue, respectively, the final effect of PTT was analyzed using the effective apoptosis retention ratio θ * eff 32 , defined as the ratio of θ * A to θ * H . Through this variable, the final PTT effect was quantitatively evaluated to identify the optimal treatment conditions. The Fig. 8 plots θ * eff as a function of φ drr and P l for various numbers of GNPs injections. The laser was set to a wavelength of 1064 nm with a gaussian profile. A similar trend to that of θ * A was observed with number of injections and it was confirmed that, for a given number of injections and φ drr , there exists a P l with optimal treatment effect. As φ drr increased, it was observed that P l must decrease for θ * eff to have a maximum value. This  www.nature.com/scientificreports/ is because an increase in φ drr is favorable for maximizing the intratumoral apoptotic temperature, but causes more thermal damage to the surrounding normal tissue. Therefore, it seems that the maximum θ * eff occurs at a lower P l than the maximum θ * A . For a φ drr of 0.2 and 0.4, the difference in θ * eff with the number of injections was negligible. This is because the distribution radius of GNPs in the tumor tissue is very small, so that an increase in the number of injections makes a negligible difference in the overall ratio. Additionally, similar to the results for θ * A , a single injection of GNPs resulted in a relatively low θ * eff compared to greater numbers of injections. This is because, in a single injection, the distribution of GNPs is concentrated in the center of the tumor tissue, hence only the temperature near the center rises to reach or exceed the apoptotic temperature range. After analyzing all cases, it was found that θ * eff is maximized when φ drr = 1, the number of injections = 7, and P l = 52 mW. These the treatment conditions produce the optimal therapeutic effect in theory, however, in the actual treatment situation, if the values of θ * eff do not differ significantly, treatment conditions that maximizes θ * A to cause more tumor tissue death may be preferable, even at the expense of thermal damage to the surrounding normal tissue, from the perspective of cancer cell death.

Conclusions
In this study, the effectiveness of PTT based on the number of GNPs injections and GNPs distribution radius in tumor tissue was quantitatively determined through numerical analysis. Numerical simulation modeling was implemented to calculate the temperature distribution inside the medium for a situation in which SCC was assumed to have developed inside a skin layer consisting of four layers and the tumor tissue, partially injected with GNPs, was irradiated with a laser. The behavior of the laser particles transmitted inside the biological tissue was analyzed using the Monte Carlo method and the optical properties of the used GNPs were calculated by DDA method.
Numerical simulations were performed to determine the temperature distribution in the medium for a given total number of GNPs while varying the number of GNPs injections into the tumor tissue, distribution radius of GNPs, and the intensity of the laser. Three metrics were used to quantitatively determine the therapeutic effect of PTT: the apoptosis retention ratio, a variable that indicates how long the temperature inside the tumor www.nature.com/scientificreports/ tissue is maintained within 43-50 °C (the temperature range over which apoptosis occurs); the thermal hazard retention value, a variable that quantitatively verifies the amount of thermal damage to the surrounding normal tissue; and the effective apoptosis retention ratio, a variable combining the former two metrics. For the given situation, the optimal treatment effect was observed when the distribution radius ratio of the injected GNPs was 1, the number of injections was 7, and the intensity of the irradiated laser was 52 mW. However, as mentioned in Section "PTT effect analysis and optimal treatment conditions", in the actual treatment, it may be necessary to choose conditions that maximize the expression of apoptosis in the tumor tissue at the expense of thermal damage to the surrounding normal tissue. In addition, since this study was only based on numerical analysis, it is considered necessary to verify it through future experimental studies. Lastly, although the temperature range at which apoptosis occurs is known, the maintenance of that temperature range is still under research and the temporal influence of apoptosis remains to be determined.